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We discuss the propagation of dark matter perturbations with non-zero velocity dispersion in cos- 
mological models. In particular a non-zero massive neutrino component may well have a significant 
effect on the matter power spectrum and cosmic microwave background anisotropy. We present a 
covariant analysis of the evolution of a dark matter distribution via a two-dimensional momentum- 
integrated hierarchy of multipole equations. This can be expanded in the velocity weight to provide 
accurate approximate equations if the matter is non-relativistic, and we also perform an expansion 
in the mass to study the propagation of relativistic matter perturbations. We suggest an approxima- 
tion to the exact hierarchy that can be used to calculate efficiently the effect of the massive neutrinos 
on the CMB power spectra. We implement the corresponding scalar mode equations numerically 
achieving a considerable reduction in computation time compared with previous approaches. 



I. INTRODUCTION 



Current observational evidence suggests that a significant fraction of the matter density of the universe is in the 
form of dark matter. Most of this must have low velocity dispersion to be consistent with large scale structure, and 
the cold dark matter (CDM) model has so far proved to be a good model. However it is quite possible that the 
dark matter has non-zero velocity dispersion, in the form of warm dark matter. There is also good evidence that 
the neutrinos have non-zero mass, so there will be a small component of hot dark matter in the form of massive 
neutrinos. In this paper we study the evolution of dark matter perturbations with a non-zero velocity dispersion, 
focusing especially on the massive neutrinos and their effect on the Cosmic Microwave Background (CMB) anisotropy. 

One of the useful products of forthcoming high-resolution CMB observations should be good limits on the massive 
neutrino masses. Atmospheric neutrino oscillation observations currently provide three-sigma evidence for a mass 
difference /S.m? ~ (3 ± 2) x 10"'^ eV^ [Q. This is consistent with a predominantly v^-Vt oscillation though more 
complicated interpretations are possible. Solar neutrino data show a much smaller mass difference Am^ (3 ± 1) x 
10"'' eV^ 1^, suggesting that two of the neutrino mass eigenstates are likely to be nearly degenerate or very small. 
The neutrino mass is related to the ratio of the density of massive neutrinos to the critical density required for a 
flat universe, by (see e.g. Ref. 0) 



.2 _ E "^i- 



n^h' = ^ " . (1) 

93.8 eV ^ ' 

For h ~ 0.7 the neutrino oscillations therefore provide the bound that 17^ > 10^"^. Recently there has also been 
controversial evidence for double beta decay [ij, which could place interesting lower limits on the neutrino mass 
eigenstates if confirmed. Current cosmological evidence is not very constraining, suggesting only that < 4.2 eV 

and hence r^^. < 0.1 §. 

If ^ 10^^ the observable effect on the cosmology is limited, with the CMB temperature and polarization power 
spectra changing by less than a percent relative to the corresponding pure CDM model. However at slightly larger 
values the effects can become important. For 57^ ^ 3 x 10^"^ there are changes in the electric polarization power 
spectrum at the few percent level. For Slj^ > 5 x 10^"^ the effect on the temperature power spectrum becomes 
significant at the percent level. Though it might be more natural for all the neutrino masses to be very small, three 
neutrinos of approximately degenerate mass (Am < 0.07 eV) would be consistent with current data and could easily 
be cosmologically interesting. 

Even if the contribution is small, it is essential to be able to include the effect of massive neutrinos in the theoretical 
calculations, both to constrain the neutrino masses themselves, and for accurate recovery of the other cosmological 
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parameters. Even if neutrino experiments suggest that the masses are very smaU we would also like to check inde- 
pendently that this is consistent with cosmology. 

In this paper we analyse the propagation of dark matter perturbations, presenting evolution equations that can be 
used for fast numerical implementations. We shall focus on linearized massive neutrino perturbations, though our 
results could be applied equally well to other forms of dark or interacting non-massless matter. In particular the 
low-energy equations could be used for the efficient propagation of warm dark matter perturbations. We also provide 
exact equations that may be a useful starting point for studying non-linear evolution. 

We use a covariant expansion of the dark matter distribution function following Ref. giving a set of exact 
multipole equations for the propagation of coUisionless matter. To compute the matter power spectrum and CMB 
anisotropies in models with hot or warm dark matter we linearize the multipole equations about an exact Friedmann- 
Robertson- Walker (FRW) model. The covariant analysis yields sets of physically transparent equations that can 
be split into scalar, vector and tensor parts as required for numerical implementation (for previous applications of 
this approach see Refs. § 0, g We derive useful limiting cases when the matter is highly relativistic and non- 
relativistic, and also discuss an accurate approximate scheme that interpolates between the two regimes. Though we 
present a covariant analysis here, our approximation techniques could be applied equally well using metric variables, 
for example the synchronous or Newtonian gauge formalism of Ref. JlOf , or the 'gauge-ready' formalism of Ref. [ pT| . 

For parameter estimation and hypothesis testing with forthcoming CMB data a large number of accurate theoretical 
power spectra will need to be generated, even if fast Monte Carlo sampling techniques are employed. This can be very 
time consuming even using efficient codes such as CMBFAST or parallelized derivatives such as CAMB ||l3|. Using 
current methods the massive neutrino evolution dominates the computation time, so our methods for speeding up the 
neutrino evolution should be useful. Our numerical code can generate CMB power spectra (including polarization) 
two or three times faster than CMBFAST on one processor, and in under a second on a modern 16 processor machine 
(non-flat models take about twice as long). 

We employ general relativity with signature (H ) and use natural units where c = 1. 



II. COVARIANT MULTIPOLE EQUATIONS 



We analyse the propagation of the dark matter with respect to a reference 4- velocity u". The variables we use are 
then physically meaningful and can in principle be measured directly by an observer moving with this velocity. We 
decompose the 4-momentum of a particle of mass m with respect to u"^ as 

= Su" + Ae°, (2) 

where E is the energy, A = V — m? is the 3-momentum and e°ea = —1. We shall assume that the dark matter is 
coUisionless so that the distribution function / = f{x°',p°') obeys the Liouville equation 

drf = (3) 

along a path described by an affine parameter r (defined so that drx"^ = P°)- 

The e" dependence of the distribution fmiction / = f{x, E, e°) can be expanded in terms of irreducible multipoles 

as 

oo 

f = Y.PA,e^'^^ + + P-be''e'' + ..., (4) 

1=0 

where the tensors Fai — F(^ai...a,) are projected (orthogonal to u") symmetric and trace-free. The angle brackets 
denote the projected symmetric trace free (PSTF) part of the enclosed indices. This multipole expansion is the 
covariant equivalent of the usual spherical harmonic expansion. Substituting into the Liouville equation we have 

oo 

{dEFAAEe^' +p^yaFA,e''^ + IFa.p'^V ae''' e^'-') = 0. (5) 

The propagation is governed by the geodesic equation p^VaP^ = 0, and the component in the u° direction gives 

{E/\)drE = drX = E^Aac" + EXiaabC^e^ - H). (6) 

This then implies that 

h\p'V,e'' = -^(A-^ + eM^e-^) - E{e''iOb'' + afeec'e^e'^ + a%e^), (7) 
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where hab is the projection tensor into the space orthogonal to Here, the local Hubble parameter, acceleration, 
shear and vorticity are given by 



Ua = U WbUa, 



D 



(0^*6), ^ab = D[QUb], 



(8) 



and T>aXa^,,,ai = hj'ha^^^ ■ ■ ■ hai^^'^bXbi...bi is the spatial covariant derivative. Using these relations in the Liouville 
equation and equating irreducible PSTF terms one obtains a set of equations for the evolution of the distribution 
function multipoles: 



+ [xEOeF^a,^^ -{I- Aa,) -^[{1 + 2)^FaA, + XEdEFaA 

I ' 2\^dEFa - "I- g I (l+l){l+2) 



A" 



21+3 [3FFa(A,_, 



(2;+3)(2;+5) [(' + ^)FFabA, + X^'^dsFabA,] cr" 

+ [X'dEFf^A,^^ -{I- 2)SF(^,_,] aa,_,a,) = 0, (9) 



where ±Fai — h^'^ Ai'u-'^'^ aFs, is the projected time derivative. 

The dark matter couples to other matter via gravity according to general relativity. The relevant quantities for the 
gravitational interaction are not the distribution function multipoles but rather the integrated quantities that make up 
the stress-energy tensor. Following Ref. |^ we therefore integrate over momentum defining the momentum-integrated 
multipoles 



~ {2l + iy. 



dXX'EFA. I ^ 



(10) 



where n = Z -I- 2i is the velocity weight. The numerical factor is introduced so that the energy density, momentum 
density, anisotropic stress, and isotropic pressure that make up the stress-energy tensor are simply 



1 7(1) 



(11) 



Integrating Eq. (|^) over momentum we obtain the set of exact multipole equations 
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This exact two-dimensional infinite hierarchy of covariant equations can be used as a starting point for analysing the 
propagation of dark matter with arbitrary velocity dispersion. In the case of coUisional matter the right hand side 
would contain a collision term. The (I, i) element of the hierarchy involves moments with / in the range (1 — 2^1 + 2) 
and i in the range (i — l,j -I- 2). The coefficients in front of the i — 1 moments vanish for i = 0, so the subset of 
equations with i > form a closed system for the i > moments. 



III. COSMOLOGICAL DARK MATTER 



The application of the two-dimensional multipole hierarchy we consider here is to dark matter perturbations in 
cosmology. The universe is assumed to be nearly FRW and we linearize by dropping products of any quantities that 
vanish in an exact FRW universe. All the I > multipoles vanish by isotropy in an exact FRW model and are 
therefore first order. The only zero-order quantities are the scale factor, Hubble parameter and the I = multipoles. 
To first order about a FRW universe the momentum-integrated multipole equations ( |l2|) become 



J 



H 



(l-n)4+i) + (3 + n)j: 



(1 + Z + n) + (3 - ; - n) [5 



I A 

3^ai 
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= 0. 



(13) 
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We characterize perturbations to the J^'^ (which do not vanish in the FRW hmit) by their first-order co-moving 
spatial gradients 

X«^5D,J«. (14) 

Here, 5* is the scale factor defined generally in a perturbed universe by integrating S = HS with initial conditions 
chosen to ensure that DaS is first order. Taking the spatial gradient of the I = equations and commuting derivatives 
gives the propagation equation 



yd) 

Aa 



(1 - 2i)J^'+^^ + {3 + 2i) J(')J ha + +h[{1~ 2^)xi'+^^ + (3 + 2i)x«J = 0, (15) 

where ha = DaS is the projected gradient of the scale factor. Note that 

ha = SiBaH - HAa) (16) 

is defined uniquely at linear order once a choice is made for u"", unlike ha itself which also depends on the choice of 
surface on which we impose S — constant as an initial condition. 

The covariant equations (|l^) and (|l^) constitute an infinite two-dimensional hierarchy for analysing the propagation 
of linear dark matter perturbations. For massless particles E — X and the velocity-weighted moments are identical, 

(i) (%') 

Jj^^ = J . The momentum-integrated equations then reduce to the usual one-dimensional Boltzmann hierarchy 
used to propagate massless neutrino perturbations Q . 

For non-relativistic matter the velocity weight I + 2i controls the magnitude of the momentum-integrated moments. 
If we neglect terms with I + 2i > the order equations become (for < 2i < n» — 2) 

+ + =0. (17) 



Using H = S/S this solves to give 



J« _ (18) 



As expected the truncated equations will be accurate to within terms that decay as an n*th-order velocity density. In 
addition to demanding that <C 1, where is the typical particle velocity, an accurate truncation also requires that 
the perturbations do not vary too rapidly in space. The spatial derivative terms in the n*th-order equations can only 
be neglected if kv^ <^ SH, i.e. the free-streaming scale in an expansion time is much less than the proper wavelength 
S/k of the fiuctuations. 

A hierarchy truncated at n» = 2 was used in Ref. fl^ to study stress effects in dark matter structure formation. 
For warm or hot dark matter one can choose as large as required to obtain accurate results. 

There is a subtlety associated with this low-velocity-weight truncation scheme that arises because the I = 1 
momentum-integrated moments are not frame-invariant in linear theory. Under linear changes in the fiducial ve- 
locity, ^ = + u", one can show that the moments transform as 



Ja! - Ja, = Ja, - [(3 + 2*) J« + (1 - 2z) (19) 

For I — 1 and i = this reduces to the well-known result q'^ = q°- — {p + p)v°'. We see that at Z = 1 the transformation 

is not homogeneous in the velocity weight. The form of Eq. (|l^) implies that the truncation condition j'^j — for 
Z -|- 2i > is only frame-invariant at linear order for n* odd. 



IV. MASSIVE NEUTRINO PERTURBATIONS 



We now discuss how to propagate massive neutrino perturbations from the early universe until the present day, 
with the aim of implementing the equations numerically in an efficient way. 

Assuming instantaneous neutrino decoupling in the early universe the zero-order distribution function for massive 
neutrinos is given by the Fermi-Dirac distribution 



(20) 
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where and Sd are the temperature and scale factor at neutrino decoupUng and q = SX is the co-moving momentum. 
We assume the particles will be highly relativistic at decoupling so the mass term can be dropped. As the universe 
expands the momenta redshift and the neutrinos evolve to become non-relativistic. 

At hydrogen recombination light neutrinos will still be quite relativistic, and the higher multipoles in the distribution 
can be significant. Usually it is assumed that the neutrinos become non-relativistic well before recombination, and 
hence the high multipoles have little effect. For light species it is necessary to evolve the higher multipoles to get 
accurate results, and approximate schemes like that in Ref. [ p^ , based on a low-Z truncation, do not give accurate 
results. For such relatively light species it is nonetheless still important to take account of the non-zero mass, as this 
can effect the matter power spectrum and CMB anisotropy power spectrum at the level of several percent (though 
some of this is trivially accounted for by the change in the background equation of state) . 

It is usual to assume that if the neutrino mass eigenstates are significant they will be approximately degenerate. 
This is reasonable as the observed mass squared differences are sufficiently small that taking them into account would 
make no observable difference to the cosmology for the foreseeable future. However it is possible that at some point 
in the future observations will be sufficiently accurate that even small mass splittings may have an observable effect, 
in which case the neutrino mass eigenstates may need to be propagated separately. 

The general regime 

When the most of the neutrinos are either highly relativistic or highly non-relativistic one can use a small mass 
or small velocity expansion, as discussed below. However in general a significant fraction of the neutrinos have 
intermediate velocities and one needs to propagate the distribution function directly to get accurate results. As the 
neutrinos are locally in equilibrium prior to decoupling, the distribution function evolves from the isotropic form in 
Eq. (pO|), but with a spatial variation in the temperature. The efficiency with which free-streaming converts this 
spatial variation into an angular variation depends on the particle velocity (and the wavelength of the fluctuation). 
For this reason the anisotropics in the distribution function do not simply inherit the momentum dependence of the 
spatial variation of the initial (isotropic) distribution, unlike the case for massless particles. 

In terms of the co-moving momentum q = SX and energy e = SE the multipole equations (^) for the propagation 
of the distribution function linearize to give 

Fa, - f |t^D^F,^, + ^B^aFA,^,) + Sn (f + if/i,,) qd^F + Sna^.a^qd^F = 0, (21) 

where the spacetime derivatives are taken at constant q. Note that 

VaFA, U = VaFA, I, + {VaSIS)qdqFA, ■ (22) 

The I — 1 equation contains the first order variable combination 

Va = S'DaF + haqdqF = S'DaF\x (23) 

that covariantly characterizes perturbations to the isotropic part of the distribution function. It integrates to give the 
co-moving gradient of the energy density: 

Xa = xi"^ ^ SBaP =-gij^ q^eVa. (24) 
Differentiating Eq. (|^) for I = with respect to time and commuting derivatives we find the propagation equation 

Va = if ^D„D''^^b + KqdqF. (25) 

This equation, together with the I > multipole equations, allows one to propagate the distribution function pertur- 
bations directly. The synchronous gauge and Newtonian gauge equivalent of these equations are well known pO| and 
are what is usually used to propagate numerically massive neutrino perturbations 0, ^]. This is slow compu- 
tationally as one has to propagate a hierarchy of multipoles for each co-moving momentum, and perform numerical 
integrations to compute the stress-energy tensor that is relevant for coupling to the other perturbations. Even for one 
massive species the numerical evolution of its distribution function dominates the computation time, and it becomes 
proportionately worse if one has to evolve several massive species. 

Note that when a given co-moving momentum becomes sufficiently non-relativistic {kq/e <C SH with k co- moving 
wavenumber) free-streaming becomes ineffective at converting spatial structure into angular structure, and the spatial- 
derivative terms in Eq. (^l]) become negligible. In this limit, FA,{q) for Z > 2 is approximately constant. When nearly 
all the particles are non-relativistic, constancy of Fa, (q) leads directly to the scaling of the momentum-integrated 
moments given in Eq. (p^). 
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The relativistic regime 

In the relativistic regime the bulk of the distribution have momenta large com par ed to the rest mass q 3> m^S. 
Defining the dimensionless mass parameter m = m^/kBT^Sd one can expand Eq. (nO) in terms of the small quantity 
fhS (which is the ratio of mass to the typical proper momentum) giving, in the FRW limit, 



where po is the density the neutrino would have if it were massless. These expressions are correct to order (fhS)'^, 
though the expansion does not generalize straightforwardly to higher order. More generally one can find a relation 
between the multipoles of different velocity weight in the perturbed universe; using the definition given in Eq. ( p^ ) 
and expanding in mS we have 



^ ^2c2 47r(-2)'(/!)^ 

^^(2^ + 1)! Jo 



jr'^j^ --^s^ oi:'::. i d,e^^^,. (27) 




Jt''^j^i^-^^S^)^J^J- (29) 



To this order we can evaluate the last integral assuming the neutrinos are massless. Starting from a locally isotropic 
distribution, the momentum dependence of the Fai follows that of the I — 1 and 2 source terms in Eq. (pi|), so that 

Fa, (X qd,F (28) 

for / > 0. Using this and integrating by parts we have 

5 

7^2 ■ 

for / > 0. In a similar manner, using Va oc qdqF in the massless limit, wc find 

X^^^ - X« (l - ^rn'S'^ « x« • (30) 

To this order in fhS, Eqs. ( |29| ) and ( ^0| ) provide closure conditions that allows reduction of the two dimensional 
momentum-integrated hierarchy, Eqs. (|l^) and (|l5|), to a one-dimensional hierarchy for propagating the J^^^ and Xa 
that are required for coupling to other perturbations. The approximation is good for (I + 2i)fh^S^ <C 1. 

The non-relativistic regime 

When the bulk of the neutrinos have q ^ m^S the energy density and pressure evolve in the FRW limit as 

180po - o , 15C5 945C7 



77r4 V 2m5 16(mS')3 ' 

_ 900po / Cs _ 63 C7 , \ 



The value of m today {S = 1) is related to the neutrino mass by 



TOi, 



1.68 X lO-'^meV. (32) 

For masses of order 0.2 eV or larger m > 10'^, and the mean speed today is less than about 10~'^c. Even in the low 
mass case the neutrinos will therefore have been non-relativistic for a significant fraction of the universe's evolution. 

Since the velocities are becoming small one can accurately truncate the momentum-integrated hierarchy, Eq. (p^), at 
some fairly low velocity weight. However in order to do this one does need to know initial values in the non-relativistic 
era for those momentum-integrated multipoles that are retained after truncation. Generally, these cannot be obtained 
from the one-dimensional hierarchy that is integrated in the relativistic regime since the eras of applicability of the 
relativistic and non-relativistic approximations do not overlap. (However, a way of interpolating between these 
two regimes with sufficient accuracy for CMB computations is discussed below.) One must therefore integrate the 
distribution function directly, using equations (|2^) and (|2^), until one enters the non-relativistic regime, at which 
point one can switch over to a truncated momentum-integrated hierarchy. 

The evolution of small-scale perturbations depends more critically on the higher velocity weight multipoles due to 
the derivative coupling in Eq. (ttq). Performing a mode expansion in wavenumber as in the appendix, the truncated 



hierarchy only becomes accurate at rather later times for small wavenumbers, as discussed in Sec. [II 
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An approximate scheme 

The relic neutrinos themselves have very low energy and will remain unobservable for the foreseeable future. We are 
therefore not interested in the neutrino perturbations per se, just in their effect on the matter and photon perturba- 
tions. In this subsection we suggest an efficient approximate scheme, based on results obtained above, for propagating 
the massive neutrino perturbations into the non-relativistic regime which we have found to be sufficiently accurate for 
the computation of the CMB power spectra at the one percent level. This approximate scheme is inadequate for the 
matter power spectrum, which is most efficiently calculated using a momentum grid in the relativistic era matched 
onto the truncated momentum-integrated hierarchy when the neutrinos are non-relativistic. 



If we use Eqs. (^6|) and (|29|), we can write the momentum-integrated hierarchy, Eqs. ( [L3D and (15), for i = during 
the relativistic era in the form 

jf^ + if [(3 + - ril - 1)] J^^^ - r^^B^,^j';^l^^ + DVfj, + SnpA,, (1 + - S^i'^i^ ~ M^a.a. = 0, (33) 

Xa + i^(r + 3)xa + D,D^jf^+3p(l + u;)/ia = 0,(34) 

where w = p/ p and r = (3w)^/^. These equations can be used to evolve the perturbations accurately while the 
majority of the particles are relativistic. 

To motivate the approximate scheme, we note that massive neutrinos perturbations affect the CMB power spectra 
in the acoustic region mainly by the damping effect on the gravitational potential of their free-streaming on these 
scales (Other effects, including changes in the position of the acoustic peaks, arise mainly from the variation in 
the background equation of state.) The effect on the potential is insensitive to the late time evolution of small-scale 
neutrino perturbations, so these can be computed in a rather coarse manner at later times. 

At late times the higher velocity weight multipoles fall off rapidly. Truncating the full momentum-integrated 
hierarchy at n* = 3, one can show that j!^^ ^ rj'i'^ and Xa^ ~ rxa^ are solutions on large scales if r « 5w. Note that 
these conditions are frame-invariant to the order of our velocity-weight truncation, and that they give the expected 
scaling as the momenta are redshifted since r oc 1/ S"^. If we use r = 5z/; in Eqs. ( pSj ) and (^4|), and drop terms with 

velocity weight n > 3, we obtain a one-dimensional hierarchy for Xa and the with / < 3 that is consistent on 
large scales with the full two-dimensional hierarchy truncated at velocity weight 3. This suggests that if we evolve 
Eqs. ( |33| ) and (^J) with r interpolating between (3z/;)^/^ and 5z/;, we obtain an approximate scheme that is accurate 
at early times and gives the correct late time behaviour for large-scale low-Z multipoles. We found that using 



mS/(mS+200) 



works well and gives CMB power spectra accurate at the one percent level, as shown in Fig. ^ The evolution of the 
neutrino modes is accurate until 3w becomes significantly less than one, as shown in Fig. ^, and is also approximately 
correct at late times on large scales. Small scales evolve incorrectly at late times, but this has no observable consequence 
for the CMB. 

The matter power spectrum is more sensitive to the late-time massive neutrino evolution than the CMB power 
spectra. On scales below the horizon size when the neutrinos become relativistic, free-streaming prevents them from 
clustering in the relativistic era and thus suppresses the growth of CDM perturbations. Once the neutrinos become 
non-relativistic they infall and cluster with the dominant CDM [l7| . The approximate scheme described above is not 
sufHciently accurate at late times for computing the matter power spectrum, and the truncated momentum-integrated 
hierarchy should be propagated instead. However the late time matter power spectrum only affects the CMB via 
second order effects such as lensing, and the CMB temperature anisotropy therefore does not depend sensitively on 
its accuracy. 



Numerical implementation 



The covariant equations can be split into scalar, vector and tensor parts and then expanded in terms of the 
appropriate eigenfunctions of the co-moving Laplacian. This gives sets of equations for the harmonic coefficients that 
can be propagated numerically. For details of the scalar and tensor cases see the appendices. 

The background evolution is affected by the neutrino density and pressure, which is a function only of fiiS. The 
numerical integrations over the distribution function required to compute background quantities therefore only need 
be performed once for each value of mS*, and the values can be stored and re-used with different neutrino masses. 
Equations (p6|) and (pll) can be used to calculate quickly the values at low and high values of fhS respectively. 
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FIG. 1: The change in the scalar CMB temperature (sohd hnes) and E polarization (dashed line) power spectra compared 
to a CDM model, with Q^h^ = 6.7 x 10"^, A^,. = 3 (top panel) and n^h^ = 2 x 10"^, iV^ = 1 (bottom panel) computed by 
integrating the distribution function perturbations. The dotted lines show the results from using the approximate Eq. (^3[), 
which agree very closely almost everywhere. We used h — 0.69, Q,a ~ 0.7, Q,k ~ 0, Q,bh^ = 0.022. 



For the computation of CMB anisotropies we require approximately one percent accuracy in the power spectrum 
for comparison with forthcoming accurate data. We implement the approximate scheme which is accurate at the one 
percent level in general, and much more so if the neutrinos are light. For one massive species this is about three times 
faster than propagating the distribution function directly, and the massive neutrino evolution no longer dominates 
the computation time. 

For comparison, and for computing the matter power spectrum, we also propagate the multipoles of the distribution 
function directly. When the momentum has redshifted sufficiently that the neutrinos are non-relativistic we switch 
to propagating the truncated momentum-integrated hierarchy. We keep terms up to = 3, giving four equations to 
propagate (the two highest order equations have analytic solutions cx S~^). The starting values for the momentum- 
integrated variables can be computed by integrating numerically over momentum at the switch-over point. Since the 
higher order variables are most important on small scales the switch-over point needs to be rather later for modes 
with large wavenumbers. This is a little slower than the approximate scheme, but can be made arbitrarily accurate 
by delaying the switch-over until the neutrinos have lower velocities. 

Hot dark matter affects the large-scale tensor power spectrum predominantly through the change in the background 
equation of state. Tensor modes interact with the neutrinos only on sub-horizon scales on which they are decaying 
away, rapidly becoming negligible compared to the scalar perturbations. 
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FIG. 2: The evolution of fc = 0.005Mpc~^ (top panel) and k = O.OlMpc"^ (bottom panel) scalar modes, with three degenerate 
massive neutrinos giving Q^h^ = 6.7 x 10~^. The blue and black dashed lines respectively show the background ionization 
fraction and 3pi,/pi,. The solid black and red lines show the massive neutrino density perturbation and momentum density, 
the dotted lines show the result of using Eq. (^^ rather than integrating the distribution function. The green and cyan lines 
show the massless neutrino and photon perturbations. The perturbations are scaled corresponding to an initial curvature 
perturbation of 0.2 and are evaluated in the zero acceleration frame. 



V. CONCLUSION 



We have derived momentum-integrated multipole equations describing the propagation of dark matter with non- 
zero velocity dispersion. In the relativistic and non-relativistic regimes there are accurate approximations that can be 
used to evolve perturbations efficiently. We focussed on massive neutrino perturbations, though the same techniques 
could be used for matter with a different distribution function. For distributions with intermediate velocities one has 
to evolve the distribution function directly to get accurate results, though we found that for computing the CMB 
power spectrum a fast approximate scheme was sufficient to take account of the effect of massive neutrinos. We have 
implemented the scalar equations numerically and our fast parallelized CMB anisotropy code is publicly available at 
http: / /camb.infc . 
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APPENDIX A: THE SCALAR EQUATIONS 



It is convenient to perform a harmonic expansion in terms of cigcnfunctions of the co-moving Laplacian 

S^Vi^ViaQ^ = fc^O^ (Al) 
We define scalar coefficients for the harmonic expansion of PSTF tensors in terms of 

QX = (5/fcyD(,, ...D,,)Q'= (A2) 

as foUows: 



(A3) 



We leave the fc-dependence of the expansion coefficients implicit. The acceleration, shear and gradient of the scale 
factor can be expanded similarly as 



Inserting the harmonic expansion into the multipole equations gives the propagation equations 



(A4) 



n 



(l-n)/f+^) + (n-3ii;(i))//*' 



i+i 

21 + 



I 7-('+l) 

21 + 1^1-1 



Si2^ka 
Sio3h' 



(3 + n)w('+i) + (l-n)w(*+2) 



SnkA 



(2 + n)w(') + (2 - 



{1 - n)w^'+^''> + (3 + n)w 



= 0, 



(A5) 



where 3w'-*' = j'*-' /p, the dash denotes the derivate with respect to conformal time, H = S' / S is the conformal 
Hubble parameter, and Pi = 1 — {P — \)K/k'^ where K is the curvature constant. 

The scalar equations for the distribution function are obtained by expanding the multipoles in terms of scalar 
harmonics as 



(A6) 



giving the multipole equations 



,^ka - Su^k (^lA 



dlnF 
ding 



0. 



(A7) 



Choosing the frame such that Aa — these equations are equivalent to the synchronous gauge equations used by 
other codes [|l^, |l2|. Assuming the neutrinos are initially highly relativistic the initial conditions for Fi{q) are the 
same as for massless neutrinos: 



//°^ dlni^ 
4 ding 



(A8) 



Later, when the particles are no longer highly relativistic, we can calculate the momentum-integrated multipoles by 
integrating over q 



' pS^ Jo 



l+2i 



FFi. 



(A9) 



We choose to use the gauge in which Aa = 0, and keep terms up to l+2i = = 3. This gives the four energy-integrated 
scalar equations 



r(0) 



' + H (4^' - 3^(1^^°)) + fc/(°' + 3 (l + h' ^ 
/W'+H(2-3«;(i))4^'+fc/(^' + 15u;(i)/i' = 



1)) + \k (2/3,4 



(0) Al) 







= 



(2 - 3«;(i)) + \k (3/334"^ - 241)) - 2fc«;(i)a = 0, 



(AlO) 
(All) 
(A12) 
(A13) 
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where 

l['Wp-'S-'. (A14) 

The density perturbation and momentum density plotted in Fig. ^ correspond to the values of Iq^'' and I^"-' in the 
Aa = frame. In the non-covariant approach Iq^^ would just be the synchronous gauge Sp/p. 

APPENDIX B: THE TENSOR EQUATIONS 

We follow the procedure in Ref. ||l^ and expand in terms of transverse tensor eigenfunctions of the co-moving 
Laplacian, 



We define scalar coefficients for the harmonic expansion of PSTF tensors with Z > 2 in terms of 



as follows: 



(Bl) 
(B2) 

(B3) 



We leave the fc-dependence of the expansion coefficients implicit. 

Inserting the harmonic expansion into the multipole equations gives the I > 2 propagation equations 



k 



(2/+l)((+l)'"' + l''i+l 2l + \'-l-\ 



ika 



(3 + n)w^'+^''> + (1 - 



0, 



(B4) 



where J^*-* = and Pi = 1 — {P — 3)K/k'^. The tensor equations for the propagation of the distribution hmction follow 
by expanding 

(2^ + 1)! 



{-2)\l\f 



giving the multipole equations 



(i+3)(i-l) a p l_„ 



2 u flln^ n 
a In q 



(B5) 



(B6) 



where Fi = 0. 
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